function [n1 n2] = func_nk_from_nz(nz, CONSTS)

    [nxy1 nxy2] = func_nxyk_from_nz(nz, CONSTS); 

    eps = CONSTS.eps;
    g = CONSTS.g;
    
    n1 = -(eps./(nz.*g)).*(nz.^2+nxy1.^2+(g^2/eps)-eps);
    
    n2 = -(eps./(nz.*g)).*(nz.^2+nxy2.^2+(g^2/eps)-eps);
    
end
